Theoretical Population Biology

Recent analyses of the effects that organismal genealogies or pedigrees of populations have on times to common ancestry for samples of genetic data are extended to cases of population subdivision and migration. Traditional coalescent models marginalize over pedigrees. A finding of a pedigree effect implies that data analysis and interpretation should not be based on the corresponding traditional coalescent model but rather on a coalescent model obtained by conditioning on the pedigree. We apply a straightforward test based on the distribution of pairwise coalescence times to four previously described scenarios of subdivision and migration. These scenarios are defined by the relative magnitudes of four parameters: the number of the local populations or demes, the deme size, the migration fraction, and the probability that migration can occur at all. We find pedigree effects in three scenarios. In two, the effect is weak if the deme size is large. The one scenario without a pedigree effect corresponds to the well known structured-coalescent model. The one scenario with a persistent pedigree effect even in the limit as the deme size tends to infinity involves long periods without gene flow interrupted by pulses of migration. We illustrate our results using simulations and numerical analysis.

Our recent work on molecular evolution and population genetics postulated that individuals with a specific mutation exhibit a fluctuation in fitness, short for FSI (fluctuating selection among individuals), whereas the fitness effect of wildtype remains a constant. An intriguing phenomenon called selection-duality emerges, that is, a slightly beneficial mutation could be a negative selection (the substitution rate less than the mutation rate). It appears that selection-duality is bounded by two bounds: the generic neutrality where the mutation is neutral by the means of fitness on average, and the substitution neutrality where the substitution rate equals to the mutation rate. In addition, the middle point of generic neutrality and substitution neutrality is called the FSI-neutrality. An important problem is about the age profile of allele frequency, i.e., the arising timing of a mutation whose frequency in the current population is given (the allele-age problem for short). Solving this problem under selection duality would help extend the standard coalescent theory that based on strict neutrality to a more general form under selection duality. In this paper, we studied the allele-age problem under selection-duality by the first arrival time approach and the mean age approach, respectively. Since the general solution of allele-age problem under selection duality is not available, we focused on solving the problem at the substitution neutrality (the up-bound of selection duality), the FSI-neutrality (the middle-point) and the generic neutrality (the low-bound), respectively. Our analysis results in an overall picture that the mean first-arrival age of a mutation at the substitution neutrality is theoretically identical to that at the FSI-neutrality, which is numerically close to that at the generic neutrality. For illustration, we calculated the mean age of nonsynonymous mutations in the human population and demonstrated that the estimated allele-age could be overestimated considerably when the effect of FSI was neglected.

The Wright-Fisher diffusion and its dual, the coalescent process, are at the core of many results and methods in population genetics. Approaches have been developed to study the dynamics of its moments under genetic drift, mutation, and recombination using ordinary differential equations. The dynamics of these moments can be used to study population genetic processes and are key building blocks of efficient methods to infer population genetic parameters, like demographic histories or fine-scale recombination rates. However, the system of equations does not close under recombination; that is, computing moments of a certain order requires knowledge of moments of higher order. By applying a coordinate transformation to the diffusion generator, we show that the canonical moments in these alternative coordinates yield a closed system, enabling more accurate numerical computations. Compared to previous approaches in the literature, we believe that this approach can be more readily extended to general scenarios. Through simulations, we verify that the derived system of differential equations can accurately capture the dynamics of the moments, and can be used to efficiently compute expected diversity and linkage statistics in population genetic samples.

Mathematically integrating genetics, development, and evolution is a longstanding challenge. Here I develop general mathematical theory that integrates sexual, discrete, multilocus genetics, development, and evolution. This yields an exact method to describe the evolutionary dynamics of allele frequencies and linkage disequilibria in multilocus systems and the associated evolutionary dynamics of mean phenotypes constructed via arbitrarily complex developmental processes. The theory shows how development affects evolution under realistic genetics, namely by shaping the fitness landscape of allele frequencies and linkage disequilibria and by constraining adaptation to an admissible evolutionary manifold (high dimensional region on the landscape) where mean phenotypes, phenotype (co-)variances, and higher moments can be developed. I derive a first-order approximation of this exact method, which yields equations in gradient form describing change in allele frequency, linkage disequilibria, and mean phenotypes as constrained, sometimes-adaptive topographies. Both the exact and approximated equations describe long-term phenotypic and genetic evolution, including the evolution of mean phenotypes, phenotype covariance matrices, "mechanistic" additive genetic cross-covariance matrices, and higher moments. I provide worked examples to illustrate the methods. The theory obtained is referred to as evo-devo dynamics, which can be interpreted as an extension of population genetics, with some similarities to quantitative genetics but with fundamental differences. The theory provides tools to re-assess empirical observations that have been paradoxical under previous theory, such as the maintenance of genetic variation, the paradox of stasis, the paradox of predictability, and the rarity of stabilising selection, which appear less paradoxical in this theory.

We study a polygenic trait under stabilizing selection at statistical equilibrium, where genetic effect, mutation rate and mutational bias are heterogeneous across loci. The model assumes L biallelic sites subject to reversible mutations, each allele described by its frequency in the population. Using a diffusion approximation, a mean-field approximation and neglecting linkage disequilibrium, we predict consistent phenomena across several regimes of selection: (1) a small deviation {Delta}* of the trait mean from its optimal value appears and persists due to genetic mutations not aligned with selection; (2) while this deviation is often undetectable at the trait level, it leaves a substantial signature at the locus level by favoring alleles reducing it, resulting in genic selection with mean coefficient s* proportional to -{Delta}* acting pervasively; (3) with stronger selection on the trait, (3a) the value of {Delta}* is decreased but the intensity of genic selection is increased in inverse proportion, resulting in an essentially constant, non negligible value of s*. We show how the stationary distribution of allelic frequencies can be obtained from {Delta}*. The latter can then be characterized as the solution to a fixed-point equation. Finally, we quantify several macroscopic observables of interest (genetic variance, description of the fluctuations of the trait mean as an Ornstein-Uhlenbeck process). The orders of magnitude of the macroscopic observables can be derived on a wide region of the parameter space. The model shows good fit and can straightforwardly be extended to accommodate pleiotropy, dominance, and some forms of epistasis. We also discuss the different breakdown which may occur (Bulmer effect, Hill-Robertson effect, breakdown of the Ornstein-Uhlenbeck approximation for the dynamics of the trait mean, depletion of genetic variability due to low mutation rates).

Estimation of the effective size (Ne) of large populations with a continuous distribution across wide geographic areas and limited dispersal of individuals has been elusive so far. Estimates of the contemporary Ne from genetic markers for such large, structured populations, typically of plant and marine species, tend to be strongly biased downwards, which has led to question their relevance. Here we show that a recently proposed estimation method of Ne from linkage disequilibrium between markers, which accounts for population structure, yields estimates of metapopulation Ne when the sampling area is sufficiently large. The method is applied to empirical data of maritime pine (Pinus pinaster Aiton). While previous estimates of Ne in pine populations were of the order of a few hundred individuals, we show that estimates of the metapopulation Ne can reach values of the order of tens of thousands of individuals. This result is especially relevant from a conservation point of view, as populations with Ne lower than 500 individuals are considered to be under the risk of extinction.

STRUCTURE and ADMIXTURE are two popular methods for detecting population structure in genetic data. They model observed genotypes as mixtures of latent ancestral populations, and the inferred admixture proportions can be used to visualize and summarize population structure. A key parameter in these models is the number of ancestral populations, K. Selecting K is a challenging problem. Perhaps the most widely used method is Evannos {Delta}K, which selects K based on the second-order change in log-likelihood as K increases. However, practitioners have often noted that {Delta}K often favors overly small K, frequently returning K = 2 even when more meaningful substructure is present. In this paper, we provide a theoretical explanation for this phenomenon: we prove that, under certain conditions, the {Delta}K method can be inconsistent, meaning that it can fail to identify the true number of populations even with infinite data.

Understanding how mutations evolve on Y chromosomes is central to explaining the origin, diversity and persistence of sex chromosomes. Mutations occurring on the Y chromosome in sexual populations experience selective dynamics that differ markedly from those on autosomes, due to a reduced effective population size and the presence of large non-recombining regions containing alleles maintained in a permanently heterozygous state. These specific features alter gene transmission in the Y chromosome population compared to autosomes, even within the same pedigree. Here, we provide a two-sex diploid Wright-Fisher model that explicitly incorporates both sex chromosomes and autosomes within a unified population framework, in order to capture the influence of these specificities on the fate of mutations, not only considering fixation probabilities but also segregation times. We use diffusion approximations and provide analytical and numerical tools to compute these quantities across a wide range of parameters and selection regimes. We recover classical results on fixation probabilities in various scenarios, including purely beneficial, deleterious or overdominant mutations, and extend them in the light of mean segregation time, a key but often overlooked determinant of evolutionary outcomes over finite timescales. In particular, our analyses show that overdominant mutations are overall more likely to fix in observable time windows on the Y chromosome than on autosomes. Individual-based simulations corroborate our approximations and highlight parameter regimes where the theoretical approach is particularly useful, especially for parameter values inducing long segregation times or small fixation probabilities, for which simulations are impractical. Our results provide a comprehensive and tractable framework for clarifying how chromosome-specific features shape evolutionary dynamics beyond fixation probabilities alone.

We study one of the simplest scenarios of polygenic selection that can be imagined: a subdivided population of diploid individuals expressing an additive trait under spatially homogeneous stabilizing selection. We are interested in the amounts of variation that can be maintained at mutation-selection-migration-drift equilibrium, at individual loci and at the level of the trait, within and among subpopulations. We derive analytical approximations for variance components and summary statistics such as FST and QST under the assumptions of the infinite-island model and compare these with individual-based simulations. We find that: (i) There is a critical migration threshold (which depends on effect sizes of trait loci) below which population structure strongly inflates genic variance in the subdivided population to levels well above those in a panmictic population. Variation within each subpopulation is maximized close to the critical migration rate. (ii) The genetic basis of trait variation across subpopulations is most similar close to this migration threshold and (counter-intuitively) decreases for higher migration rates. This has consequences for the portability of Genome-Wide Association Studies (GWAS) between subpopulations, i.e, the extent to which loci with large contributions to variance in one subpopulation explain variance in other subpopulations. (iii) An analytical mean-field approach based on the single-locus diffusion approximation, together with effective migration and selection parameters (to account for associations between loci), very accurately predicts various quantities.

We model the introgression of a genome with many weakly selected linked loci into a large homogeneous population, under the simple assumption that a block of introduced genome has a selective effect proportional to its map length. Using a diffusion approximation, we compute the probability that some part of the initial genome survives the initial phase of the introgression and give the typical length of the surviving blocks. Our results quantify the effect of recombination on selection and drift during an introgression and indicate that the fate of the genome depends on the strength of selection relative to recombination. When selection is positive some parts of the genome are able to survive at large times, but large blocks can only persist if selection is stronger than recombination. Surprisingly, the probability of such a successful introgression is independent of the strength of recombination and is the same as that for a single beneficial allele. Conversely, a deleterious or neutral genome is eventually lost, but at a much slower rate than a single allele with the same selective effect. In this case, surviving blocks are very small. We also consider the introgression of a genome made of a single beneficial allele linked to a deleterious background and compute the amount of deleterious material that hitchhikes during fixation.

Alleles with opposing effects on fitness characters are said to exhibit selectional antagonistic pleiotropy (broadly construed so that effects are not necessarily confined to the same individual). A number of theoretical investigations considered the case where a pair of alleles at a locus influences two fitness components and derived the conditions giving rise to stable polymorphism under various assumptions about the mode of trait-interaction. Strikingly, many of these analyses concluded that the potential for maintaining polymorphism is strongly constrained by the joint influence of two factors: (1) the prevalence of weak selection coefficients over coefficients of large magnitude, and (2) the absence of beneficial dominance reversals (where the deleterious effects of each allele are partially or completely masked in the heterozygous genotype). Consequently, the conclusion that selective polymorphism is unlikely to be maintained by intralocus mechanisms of antagonistic pleiotropy has achieved widespread acceptance. Here we argue that such conclusions do not apply to any of the following models of antagonism: (i) additive trait-interaction, (ii) multiplicative trait-interaction, (iii) bivoltine selection, (iv) soft selection, (v) hard selection, and (vi) sexual antagonism. We demonstrate that the parameter space giving rise to stable allelic variation is quite large throughout, and moreover, the plenitude of suitable parameters neither depends on the strength of selection nor requires dominance reversal. Dominance coefficients associated with stringent conditions for stable polymorphism are shown to be atypical as compared to all feasible parameters, and best regarded as an outcome of adherence to a special relation: dominance with a constant magnitude and direction, which includes the case of additive allelic effects at a locus. Properties of single-locus equilibria (heterozygosity, allele frequency differentiation) are investigated, as well as the contribution of dominance schemes to the genetic variance in fitness characters in populations at multilocus linkage equilibrium. Author summaryAllelic variants at a locus with opposing effects on multiple fitness components (antagonistic fitness pleiotropy) have long been appreciated as a possible source of balancing selection. The prevalence of polymorphism owing to this form of natural selection, however, has been doubted on theoretical grounds due to the fact that standard assumptions of genetic models (namely, constant magnitudes for the dominance coefficients) are hardly conducive to the maintenance of polymorphism. The major exception to this conclusion lies with schemes that exhibit dominance reversal (where the direction of dominance for antagonistic alleles flips across fitness components). Here we conduct a geometric analysis of the space of polymorphism-promoting dominance parameters and conclude that the conditions for maintaining balanced alleles is unrestrictive, with non-reversals playing an underappreciated role.

The spatial structure of populations may promote the emergence and maintenance of cooperation. Cooperation in the prisoners dilemma is favored under specific update rules in evolutionary graph theory models with one individual per node of a graph, but this effect vanishes in models with well-mixed demes connected by migrations under soft selection. In contrast, experiments and models involving cycles of growth, merging and dilution have shown that spatial structure can favor cooperation. Here, we reconcile these findings by studying deme-structured populations under growth-merging-dilution dynamics, corresponding to a clique (fully connected graph) under hard selection. We obtain analytical conditions for the cooperator fraction to increase during deterministic logistic growth, and to increase on average under dilution-growth-merging cycles, in the weak selection regime. Furthermore, we analytically express the fixation probability of cooperators under weak selection, yielding a criterion for cooperative mutants to have a higher fixation probability than neutral ones. Finally, numerical simulations show that stochastic growth further promotes cooperation. Overall, hard selection is essential for cooperation to be promoted in deme-structured populations.

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

A common assumption in molecular evolution is the fixed selection nature of a mutation, for instance, a neutral mutation is selectively neutral for all individuals who carry the mutation, and so forth a deleterious or beneficial mutation. Our recent work challenged this presumption, postulating that individuals with a specific mutation exhibit a fluctuation in fitness, short for FSI (fluctuating selection among individuals). Moreover, an intriguing phenomenon called selection-duality emerges, that is, a slightly beneficial mutation could be a negative selection (the substitution rate less than the mutation rate). It appears that selection-duality is bounded: the low-bound is the generic neutrality where the mutation is neutral by the means of fitness on average, while the up-bound is the substitution neutrality where the substitution rate equals to the mutation rate. In this paper, we conducted a thorough theoretical analysis to evaluate how many generations needed for a selection-duality mutation to be fixed in a finite population. A striking finding is that the mean fixation time of a selection-duality mutant, including the generic neutrality and the substitution neutrality, is approximately identical, which is considerably shorter than the case of strict neutrality without FSI. One may further envisage that the fast-fixation nature of selection-duality mutations could result in a considerable genetic reduction at linked sites.

The analysis of large gene and metabolic networks is often hindered by unknown biochemical parameters and the nonlinear nature of classical S-system models. To address this, we introduce a framework based on combinatorial toric geometry computed with tools such as Normaliz, SageMath, it is worth mentioning this technique in not restrictive to integer vectors, there exists a natural extension to real geometries. Unlike traditional approaches, which rely on parameter dependent fixed points, our method constructs a Topological Environment derived from the dual space of kinetic orders, leading to what we call orthogonal enzyme kinetics. Within this topological setting, fixed points are computed on the algebraic torus, enabling the transformation of nonlinear dynamics into linear forms. Importantly, these fixed points are independent of kinetic parameters and depend only on network topology and interaction signs. Applying this methodology to gene circuits involved in circadian rhythms, we reproduce previously reported oscillatory physiologies.

Males during mating never transfer just sperm; to the best of our knowledge, they always deliver a rich seminal fluid as well. The proteins in the seminal fluid have an important sperm supporting role, but they also cause changes in the female physiology and can impose a mating cost. The associated costs and delay in the time of remating, lead to the view that those proteins evolved primarily due to sexual conflict and delay female remating beyond the optimal rate. To examine the role of seminal fluid proteins in sexual conflict we use a mathematical model of reproductive physiology, informed by the accumulated knowledge on Drosophila melanogaster. In accordance with the theory, we find that males always benefit from inducing longer remating intervals in females. But, we also find that this conflict is reduced when female reproduction is regulated by the male proteins. Without seminal fluid proteins females have a single, well-defined, optimal remating rate. However, when seminal proteins are used to regulate females reproduction, females can reach the same offspring production for a range of mating intervals. This wider range of possible remating times could provide females with a buffer against uncertain mating opportunity. It could also allow females to be more selective on male quality, by reducing the cost associated with delaying remating. Our results suggest that, while there is a conflict over the remating rate, seminal fluid proteins reduce its intensity, highlighting their role in aligning the interests of both sexes.

Gene tree parsimony (GTP) is a common approach for efficient reconciliation of multiple discordant gene tree phylogenies for the inference of a single species tree. However, despite the popularity of GTP methods due to their low computational costs, prior work has shown that some commonly employed parsimony costs are statistically inconsistent under the multispecies coalescent process. Furthermore, a fine-grained analysis of the inconsistency has indicated potentially complimentary behavior of duplication and deep coalescence costs for symmetric and asymmetric species trees. In this work, we prove inconsistency of GTP estimators for all linear combinations of duplication, loss and deep coalescence scores. We also explore empirical implications of this result evaluating inference results of several GTP cost schemes under varying levels of incomplete lineage sorting.

Bayesian inference of divergence times for extant species using molecular data is an unconventional statistical problem: Divergence times and molecular rates are confounded, and only their product, the molecular branch length, is statistically identifiable. This means we must use priors on times and rates to break the identifiability problem. As a consequence, there is a lower bound in the uncertainty that can be attained under infinite data for estimates of evolutionary timescales using the molecular clock. With infinite data (i.e., an infinite number of sites and loci in the alignment) uncertainty in ages of nodes in phylogenies increases proportionally with their mean age, such that older nodes have higher uncertainty than younger nodes. On the other hand, if extinct taxa are present in the phylogeny, and if their sampling times are known (i.e., heterochronous data), then times and rates are identifiable and uncertainties of inferred times and rates go to zero with infinite data. However, in real heterochronous datasets (such as viruses and bacteria), alignments tend to be small and how much uncertainty is present and how it can be reduced as a function of data size are questions that have not been explored. This is clearly important for our understanding of the tempo and mode of microbial evolution using the molecular clock. Here we conducted extensive simulation experiments and analyses of empirical data to develop the infinite-sites theory for heterochronous data. Contrary to expectations, we find that uncertainty in ages of internal nodes scales positively with the distance to their closest tip with known age (i.e., calibration age), not their absolute age. Our results also demonstrate that estimation uncertainty decreases with calibration age more slowly in data sets with more, rather than fewer site patterns, although overall uncertainty is lower in the former. Our statistical framework establishes the minimum uncertainty that can be attained with perfect calibrations and sequence data that are effectively infinitely informative. Finally, we discuss the implications for viral sequence data sets. In a vast majority of cases viral data from outbreaks is not sufficiently informative to display infinite-sites behaviour and thus all estimates of evolutionary timescales will be associated with a degree of uncertainty that will depend on the size of the data set, its information content, and the complexity of the model. We anticipate that our framework is useful to determine such theoretical limits in empirical analyses of microbial outbreaks.

Organisms must manage a trade-off between robustness and flexibility as they enact adaptive behaviors. One way organisms achieve this is by navigating a network of quasi-stable behavioral states. Evidence for such behavioral states has been observed in many organisms, and new methods for detecting these states have taken on a prominent research focus. Although dynamical models demonstrating behavioral switching have been developed significantly over the past few decades, theories of the similarities and differences among these models, necessary for advancing empirical modeling, have not yet been fully elaborated. Here, we consider behavioral switching in two different classes of dynamical models of the forward-reversal behavioral transition in C. elegans. We first show how fundamentally different models can give rise to similar phenomena under noisy conditions. We then analyze the deterministic aspects of these models to expand on their differences, clarifying the theoretical relationship between them. Finally, we demonstrate how sequence models can be further extended to incorporate dwell times for behavioral states. Our work contributes toward a broader theoretical understanding of behavioral switching in adaptive systems.

Article summaryGene interactions are common, yet additive genetic models often predict short-term evolution and breeding response. This study argues that additivity can arise because populations sample only a small neighbourhood of a curved fitness landscape. In additive channels, genetic variation is small enough that local curvature contributes little to heritable fitness differences. The study defines an additivity index ([A]g) that compares variance from the local slope of log-fitness with variance from curvature, and links this ratio to expected prediction accuracy under Gaussian assumptions. A selection-inheritance framework shows when additive channels persist and when populations leave them. It yields testable predictions.